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Abstract 

We introduce a finite element construction for use on the class of convex, planar polygons and show 
it obtains a quadratic error convergence estimate. On a convex n-gon satisfying simple geometric crite- 
ria, our construction produces 2n basis functions, associated in a Lagrange-like fashion to each vertex 
and each edge midpoint, by transforming and combining a set of n(n + l)/2 basis functions known to 
obtain quadratic convergence. The technique broadens the scope of the so-called 'serendipity' elements, 
previously studied only for quadrilateral and regular hexahedral meshes, by employing the theory of 
generalized barycentric coordinates. Uniform a priori error estimates are established over the class of 
convex quadrilaterals with bounded aspect ratio as well as over the class of generic convex planar polygons 
satisfying additional shape regularity conditions to exclude large interior angles and short edges. Nu- 
merical evidence is provided on a trapezoidal quadrilateral mesh, previously not amenable to serendipity 
constructions, and applications to adaptive meshing are discussed. 

1 Introduction 

Barycentric coordinates provide a basis for linear finite elements on simplices, and generalized barycentric 
coordinates naturally produce a suitable basis for linear finite elements on general polygons. Various appli- 
cations make use of this technique [15, 16, 25, 26, 28, 30, 32, 33, 34, 37], but in each case, only linear error 
estimates can be asserted. A quadratic finite element can easily be constructed by taking pairwise products 
of the basis functions from the linear element, yet this approach has not been pursued, primarily since the 
requisite number of basis functions grows quadratically in the number of vertices of the polygon. Still, many 
of the pairwise products are zero along the entire polygonal boundary and thus are unimportant for inter- 
element continuity, a key ingredient in finite element theory. For quadrilateral elements, these 'extra' basis 
functions arc well understood and, for quadrilaterals that can be affinely mapped to a square, the so-called 
'serendipity element' yields an acceptable basis consisting of only those basis functions needed to guaran- 
tee inter-element continuity [39, 4, 3]. We generalize this construction to produce a quadratic serendipity 
element for arbitrary convex polygons derived from generalized barycentric coordinates. 

Our construction yields a set of Lagrange-like basis functions {ipij} - one per vertex and one per edge 
midpoint - using a linear combination of pairwise products of generalized barycentric functions {Aj}. We 
show that this set spans all constant, linear, and quadratic polynomials, making it suitable for finite clement 
analysis via the Bramble-Hilbert lemma. Further, given uniform bounds on the aspect ratio, edge length, 
and interior angles of the polygon, we bound llV'ijlliji(Q) uniformly with respect to jlAjjlji^i^fj-). Since our 
previous work shows that jjAill^i^ji^ is bounded uniformly under these geometric hypotheses, this proves 
that the ipij functions arc well-behaved. 
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Figure 1: Overview of the construction process. In each figure, the dots are in one-to-one correspondence 

with the set of functions listed below it. At filled dots, all functions in the set evaluate to zero except for the 
function corresponding to the dot which evaluates to one. The rightmost element has quadratic precision 
with only these types of 'Lagrange-like' basis functions. 



Figure 1 gives a visual depiction of the construction process. Starting with one generalized barycentric 

fimction per vertex of an n-gon, take all pairwise products yielding a total of n(n + l)/2 functions 
Mafe ■= AflAb. The linear transformation A reduces the set {j^ab} to the 2n element set indexed over 

vertices and edge midpoints of the polygon. A simple bounded linear transformation B converts {^ij} into a 
basis {ilJij} which satisfies the "Lagrange property" meaning each function takes the value 1 at its associated 
node and at all other nodes. 

The paper is organized as follows. In Section 2 we review relevant background on finite element theory, 
serendipity elements, and generalized barycentric functions. In Section 3, we show that if the entries of 
matrix A satisfy certain linear constraints Q^l-Q^3, the resulting set of functions {^ij} span all constant, 
linear and quadratic monomials in two variables, a requirement for quadratic finite elements. In Section 4, we 
show how the constraints Q^l-Q^S can be satisfied in the special cases of the imit square, regular polygons, 
and convex quadrilaterals. In Section 5, we show how Q^l-Q^S can be satisfied on a simple convex polygon. 
We also prove that the resulting value of ||A|| is bounded uniformly, provided the convex polygon satisfies 
certain geometric quality conditions. In Section 6 we define B and show that the final {ipij} basis is Lagrange- 
like. Finally, in Section 7, we describe practical applications, give numerical evidence, and consider future 
directions. 




Figure 2: Notation used to describe polygonal geometry. 
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2 Background and Notation 



Let Q be a convex polygon with n vertices (vi, . . . , v„) ordered counter-clockwise. Denote the interior angle 
at Vi by f3i. The largest distance between two points in fl (the diameter of fl) is denoted diam(J7) and the 
radius of the largest inscribed circle is denoted p{il). The center of this circle is denoted c and is selected 
arbitrarily when no unique circle exists. The aspect ratio (or chunkiness parameter) 7 is the ratio of the 
diameter to the radius of the largest inscribed circle, i.e. 

diani(r2) 

The notation is shown in Figure 2. 

For a multi-index a = (ai,a2) and point x = {x,y), define x" := x"^y"^, a\ := aia2, \a\ := ai + a2, 
and D"u := d^'^^u/dx'^^dy'^K The Sobolev semi-norms and norms over an open set J7 are defined by 

|w|H-(n) := / ^ \D°'u{yi)\^ dyi and I|w||h-(q) — |w|^'=(Q) • 

\a\=m a<k<m 

The ff°-norm is the L-^-norm and will be denoted ||- 11^2(0). The space of polynomials of degree < fc on a 
domain is denoted Vk- 



2.1 The Bramble-Hilbert Lemma 

A finite element method approximates a function u from an infinite-dimensional functional space by a 
function Uh from a finite-dimensional subspace Vh <Z V. One goal of such approaches is to prove that the 
error of the numerical solution Uh is bounded a priori by the error of the best approximation available in 
Vh, i.e. ||u — UhWy < Ciniyj^Vh 11^ ~ '"^lly this paper, V = and Vh is the span of a set of functions 
defined piecewise over a 2D mesh of convex polygons. The parameter h indicates the maximum diameter 
of an element in the mesh. Further details on the finite element method can be found in a number of 
textbooks [8, 5, 11, 39]. 

A quadratic finite element method in this context means that when /i ^ 0, the best approximation error 
(inf^gVh II" ^ "'IIv) converges to zero with order h?. This means the space Vh is 'dense enough' in V to 
allow for quadratic convergence. Such arguments are usually proved via the Bramble-Hilbert lemma which 
guarantees that if Vh contains polynomials up to a certain degree, a bound on the approximation error can be 
found. The variant of the Bramble-Hilbert lemma stated below includes a uniform constant over all convex 
domains which is a necessary detail in the context of general polygonal elements and generalized barycentric 
functions. 

Lemma 2.1 (Bramble-Hilbert [35, 10]). There exists a uniform constant Cbh such that for all con- 
vex polygons and for all u G -ff'^(fi), there exists a degree k polynomial p„ with \\u — PuWh^' (q) < 

Csi/ diam(f7)''+i~'' |u|^fc+i(Q) for any k' < k. 

Our focus is on quadratic elements (i.e., fc = 2) and error estimates in the i/^-norm (i.e., fc' ~ 1) which 
yields an estimate that scales with diam(il)^. Our methods extend to more general Sobolev spaces (i.e., 
W'^'^, the space of functions with all derivatives of order < fc in L^') whenever the Bramble-Hilbert lemma 
holds. Extensions to higher order elements (fc > 2) will be briefly discussed in Section 7. 

Observe that if fl is transformed by any invertible affine map T, the polynomial p o ^ on Tfl has the 
same degree as the polynomial p on Q,. This fact is often exploited in the simpler and well-studied case 
of triangTilar meshes; an estimate on a reference triangle K becomes an estimate on any physical triangle 
K by passing through an affine transformation taking K to K. For n > 3, however, two generic n-gons 
may differ by a non-affine transformation and thus, as we will see in the next section, the use of a single 
reference element can become overly restrictive on element geometry. In our arguments, we instead analyze 
classes of "reference" elements, namely, diameter one convex quadrilaterals or convex polygons of diameter 
one satisfying the geometric criteria given in Section 2.3; see Figure 3. Using a class of reference elements 
allows us to establish uniform error estimates over all afiine transformations of this class. 
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Figure 3: Using affine transformation, analysis can be restricted to a class of unit diameter polygons. 

2.2 Serendipity Quadratic Elements 

The term 'serendipity element' refers to a long-standing observation in the finite element community that 
tensor product bases of polynomials on rectangular meshes of quadrilaterals in 2D or cubes in 3D can obtain 
higher order convergence rates with fewer than the 'expected' number of basis functions resulting from tensor 
products. This phenomenon is discussed in many finite element textbooks, e.g. [31, 20, 8], and was recently 
characterized precisely by Arnold and Awanou [3]. For instance, the degree r tensor product basis on a 
square reference element has (r + 1)^ basis functions and can have guaranteed convergence rates of order 
r + 1 when transformed to a rectangular mesh via bilinear isomorphisms [4] . By the Bramble-Hilbert lemma, 
however, the function space spanned by this basis may be unnecessarily large as the dimension of Vr is only 
(r + l)(r + 2)/2 and only 4r degrees of freedom associated to the boundary are needed to ensure sufficient 
inter-element continuity in H^. 

This motivates the construction of the serendipity element for quadrilaterals. By a judicious choice of 
basis functions, an order r convergence rate can be obtained with one basis function associated to each 
vertex, (r — 1) basis functions associated to each edge, and q additional functions associated to interior 
points of the quadrilateral, where g = for r < 4 and q= (r — 2)(r — l)/2 for r > 4 [3]. Such an approach 
only works if the reference element is mapped via an afiine transformation; it has been demonstrated that 
the serendipity element fails on trapezoidal elements, such as those shown in Figure 10 [24, 22, 39, 38]. 

Some very specific serendipity elements have been constructed for quadrilaterals and regular hexagons 
based on the Wachspress coordinates (discussed in the next sections) [36, 2, 18, 1, 19]. Our work generalizes 
this construction to arbitrary polygons without dependence on the type of generalized barycentric coordinate 
selected and with uniform bounds under certain geometric criteria. 

2.3 Generalized Barycentric Elements 

To avoid non-aflane transformations associated with tensor products constructions on a single reference 
element, we use generalized barycentric coordinates to define our basis functions. These coordinates are any 
functions satisfying the following agreed-upon definition in the literature. 

Definition 2.2. Functions Aj : ^ R, i = 1, . . . ,n are barycentric coordinates on Q. if they satisfy two 
properties. 

Bl. Non-negative: Aj > on Q.. 

n 

B2. Linear Completeness: For any linear function i : O ^ R, L = L(vj)Ai. 

i=l 

We will further restrict our attention to barycentric coordinates satisfying the following invariance prop- 
erty. Let T : R^ — >■ be a composition of translation, rotation, and uniform scaling transformations and 
let {Af } denote a set of barycentric coordinates on Tf2. 
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B3. Invariance: Ai(x) = Xl{T{x)). 



This assumption will allow estimates over the class of convex sets with diameter one to be immediately 

extended to generic sizes since translation, rotation and uniform scaling operations can be easily passed 
through Sobolev norms. At the expense of requiring uniform bounds over a class of diameter-one domains 
rather than a single reference element, we avoid having to handle non-afiine mappings between reference and 
physical elements. 

A set of barycentric coordinates {Aj} also satisfies three additional familiar properties. A proof that Bl 
and B2 imply the additional properties B4-B6 can be found in [17]. Note that B4 and B5 follow immediately 
by setting L = 1 or L = x in B2. 

n 

B4. Peirtition of unity: Aj = 1. 

i=l 

n 

B5. Linear precision: ^^ViAi(x) = x. 

i=l 

B6. Interpolation: Ai(vj) = Sij. 

Various particular barycentric coordinates have been constructed in the literature. We briefly mention a 
few of the more prominent kinds and associated references here; readers arc referred to our prior work [17, 
Section 2] as well as the survey papers of Cueto et al. [9] and Sukumar and Tabarraei [33] for further 
details. The triangulation coordinates A"^' are defined by triangulating the polygon and using the standard 
barycentric coordinates over each triangle [14]. Harmonic coordinates A^'*'' arc defined as the solution to 
Laplace's equation on the polygon with piecewise linear boundary data satisfying B6 [21, 25, 7]. Explicitly 
constructed functions include the rational Wachspress coordinates A^*'''^ [36], the Sibson coordinates A^'*^^ 
defined in terms of the Voronoi diagram of the vertices of the polygon [29, 12], and the mean value coordinates 
^MVai defined by Floater [13, 14]. 

To obtain convergence estimates with any of these functions, certain geometric conditions must be satisfied 
by a generic mesh element. We will consider domains satisfying the following three geometric conditions. 

Gl. Bounded aspect ratio: There exists 7* e M such that 7 < 7*. 

G2. Minimum edge length: There exists rf* e M such that [vj — Vj| > rf, > for all i ^ j. 
G3. Maximum interior angle: There exists /3* e M such that Pi < j3* <-k for all i. 

Under some set of these conditions, the H^-noxvn of many generalized barycentric coordinates are bounded 
in norm. This is a key estimate in asserting the expected (linear) convergence rate in the typical finite 
element setting. 

Theorem 2.3 ([27] for A^^^ai j^yj f^^. others). For any convex polygon O satisfying Gl, G2, and G3, 
A™, A^ar, j^w^ach^ j^sibs^ ;s^MVai i^g^rided m H^, i.e. there exists a constant C > such that 

||Ai||jji(Qj < C. (2.1) 

The results in [17] and [27] are somewhat stronger than the statement of Theorem 2.3, namely, not all of 
the geometric hypotheses are necessary for every coordinate type. Our results, however, rely generically on 
any set of barycentric coordinates satisfying (2.1). Any additional dependence on the shape geometry will 
be made explicitly clear in the proofs. Weakening of the geometric hypotheses is discussed in Section 7. 
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2.4 Quadratic Precision Barycentric Functions 

Since generalized barycentric; coordinates are only guaranteed to have linear precision (property B5), they 
cannot provide greater than linear order error estimates. Pairwise products of barycentric coordinates, 
however, provide quadratic precision as the following simple proposition explains. 

Proposition 2.4. Given a set of barycentric coordinates {A,}-Li, the set of functions {nab} '■= {^a^b}2 b=i 
has constant, linear, and quadratic precision^, i.e. 

n n n n n n 

^^l^ab = '^, ^^^alJ'ab = '>^ and ^^VaV^/i„6 = XX^. (2.2) 

a=l 6=1 a=l 6=1 a=l 6=1 

Proof. The result is immediate from properties B4 and B5 of the Aj functions. □ 

The product rule ensures that Theorem 2.3 extends immediately to the pairwise product functions. 

Corollary 2.5. Let Q be a convex polygon satisfying Gl, G2, and G3, and let Xi denote a set of barycentric 
coordinates satisfying the result of Theorem 2.3 (e.g. A^^'', A^^'^'^, A^'^^, or X^^'^^). Then pairwise 
products of the Aj functions are all bounded in H^, i.e. there exists a constant C > such that 

llMa6||ffi(n) ^ C'- (2-3) 

While the {/iafe} functions are commonly used on triangles to provide a quadratic Lagrange element, 
they have not been considered in the context of generalized barycentric coordinates on convex polygons as 
considered here. Langer and Seidel have considered higher order barycentric interpolation in the computer 
graphics literature [23]; their approach, however, is for problems requiring C^-continuous interpolation rather 
than the weaker _ff^-contimiity required for finite element theory. 

In the remainder of this section, we describe notation that will be used to index functions throughout 
the rest of the paper. Since /lab = l-^ba, the summations from (2.2) can be written in a symmetric expansion. 
Define the paired index set 

/:= {{a, 6} I a, 6 e {!,..., n}}. 

Note that sets with cardinality 1 occur when a = b and are included in /. We partition / into three subsets 
corresponding to geometrical features of the polygon: vertices, edges of the boundary, and interior diagonals. 
More precisely, I = V l^ El^ D, a disjoint union, where 

V := {{a,a}\ae {!,..., n}} ; 
E := {{a,a + l}\ae {!,... ,n}}] 
D := I\{VUE). 

In the definition of E above (and in general for indices throughout the paper) , values are interpreted modulo 
n, i.e. {n,n + 1}, {n, 1}, and {0, 1} all correspond to the edge between vertex 1 and vertex n. To simplify 
notation, we will omit the braces and commas when referring to elements of the index set /. For instance, 
instead of M{a,6} , we write just Hab- We emphasize that ab & I refers to an unordered and possibly non-distinct 
pair of vertices. Occasionally we will also use the abbreviated notation 

Va + Vft 
Va6 := ^ , 

SO that Vaa is just a different expression for Vq. Under these conventions, the precision properties from (2.2) 
can be rewritten as follows. 



^Note that xx^ is a symmetric matrix of quadratic monomials. 
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Ql. Constcint Precision: /i^a + ^ 2/ia6 = 1 

aaGV abeEUD 

Q2. Linear Precision: ^ ^ '^aai^-aa + 

aaeV abe-EUD 

Q3. Quadratic Precision: ^ ^a^llJ-aa + ^ (vaV^ + Vbvl)i^ab = xx"^ 

aaey abe-EUD 

3 Reducing Quadratic Elements to Serendipity Elements 

We now seek to reduce the set of pairwise product functions {Hab} to a basis {^ij} for a serendipity quadratic 
finite element space. Our desired basis must 

(i) span all quadratic polynomials of two variables on 

(ii) be exactly the space of quadratic polynomials (of one variable) when restricted to edges of il, and 

(iii) contain only 2n basis functions. 

The intuition for how to achieve this is seen from the number of distinct pairwise products: 
\M\ = \I\ = \V\ + \E\ + \D\=n + n+ ^^^^ = " + (2) 

On dn, functions with indices in V vanish on all but two adjacent edges, functions with indices in E vanish 
on all but one edge, and functions with indices in D vanish on all edges. Since Q1-Q3 hold on all of J7, 
including dil, the set {/j-ab : ab G V U E} satisfies (ii) and (iii), but not necessarily (i). Thus, our goal is 
to add linear combinations of functions with indices in D to those with indices in V or E such that (i) is 
ensured. 

We formalize this goal as a linear algebra problem: find a matrix A for the equation 

:= A[^Iab] (3.1) 
such that l^ij] satisfies the following conditions analogous to Q1-Q3: 

Qjl. Constant Precision: S^u + 2^j(j_|_i) = 1. 

iiev ,(i+i)e£ 

Q|2. Linear Precision: ^ vu^u + ^ 2vj(j+i)^j(j+i) = x. 

iiev 

Q|3. Quadratic Precision: 

X] '^i'^rCii + (^'^f+l + Vi+ivf )$i(i+i) = XX^. 

iiev 

Since (3.1) is a linear relationship, we are still able to restrict our analysis to a reference set of unit diameter 

polygons (recall Figure 3). Specifically if matrix A yields a "reference" basis T[^jj] = AT[i^ab\ satisfying 
Qjl-Q^3, then the "physical" basis [^ij] = Aluab] also satisfies Q^l-Qj3. 
To specify A in (3.1), we will use the specific basis orderings 

Kij] [ $11) $22j ■ • ■ ) Cranj $12 ) $23 ) • • • ) $(n- l)n ; $n(n+l) 

], (3.2) 

' V ' V ^ ' 

indices in V indices in E 
[Maft] := [ Mll)M22, • • • /'l2-/'2:! /' („ - 1) „ ■ /' (ri+l) , (3.3) 

^ V ' '•^ ^ ' 

indices in V indices in E 

/ii3, . . . , (lexicographical), . . . , M(n-2)n ]• 

^ V " 

indices in D 
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The entries of A are denoted c^^, following the orderings given in (3.2)-(3.3) so that 













'^(n-2)Ti 


n(n+l) 


n(n+l) 


n(n+l) 
^(n-2)n . 



(3.4) 



A sufficient set of constraints on the coefficients of A to ensure Q^l-Q^3 is given by the following lemma. 

Lemma 3.1. The constraints Q^l-Q^3 listed below imply Q^l-Q^3, respectively. That is, Q^l Q^l, 
Q^2 ^ Qj2, and Q^3 ^ Q^3. 

E C + E 24(i+'^ = 1 Vaa e V, and 
EC+ E 2cl^^)=2, VafeeEUZ?. 

iieV i{i+l)eE 

QJ- E + E 2c;(^+^^Wi(i+i) = Vaa Vaa € V, and 

iieV i(i+l)eE 

E ^"b^" + E 2cir'^ = 2t;a6, Va6 e U D. 

iieV 



E C^i^f + E 4^i+'H«i«f+i + ) = «a«r Va e V, and 

iieV i{i+l)eE 

E c:,v,vl + E 4r'Vi^r+i + Vi+ivJ) = Vavl + Vbvl, Va6 € U D. 

i(i+l)eE 



iieV 

Proof. Suppose Q^l holds. Substituting the expressions from Q^l into the coefficients of Ql (from the end 
of Section 2), we get 

E (E^aa+ E 2<(r^MMaa + 

E fE-^+ E K'rUi^ab = i. 

abeEUD \iieV i{i+l)eE 

Regrouping this summation over ij indices instead of ab indices, we have 



EfE^M+ E 2fE^:rVab)=i. 

iieV Vabe/ / ,(i+l)gE \abel ) 



(3.5) 



Since (3.1) defines — Cabl^ab, (3-5) is exactly the statement of Q^l. The other two cases follow by the 



abel 



same technique of regrouping summations. 



□ 



Wc now give sonic remarks about our approach to finding coefficients satisfying Q^l-Q^S. Observe that 
the first equation in each of Q^l-Q^S is satisfied by 



''aa • "«a "^^"-^ '^aa ■ '-' 



(3.6) 
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a{a+l 



Figure 4: When constructing the matrix A, only six non-zero elements are used in each column corresponding 

to an interior diagonal of the pairwisc product basis. In the serendipity basis, the interior diagonal function 
/Xab only contributes to six basis functions as shown, corresponding to the vertices of the diagonal's endpoints 
and the midpoints of adjacent boundary edges. 

Further, if ab = a{a + 1) G E, the second equation in each of Q^l-Q^S is satisfied by 

ciV+i):=0 and c^^^^^ := 5,, (3.7) 

The choices in (3.6) and (3.7) give A the simple structure 

A := [ I I A' ] , (3.8) 

where I is the 2n x 2n identity matrix. Note that this corresponds exactly to our intuitive approach of 
setting each function to be the corresponding function plus a linear combination of Hab functions with 
ab e D. Also, with this selection, we can verify that many of the conditions which are part of Q^l, Q^2 and 
Qg3 hold. Specifically, whenever ab G V U E, the corresponding conditions hold as stated in the following 
lemma. 

Lemma 3.2. The first 2n columns of the matrix A given by (3.8), i.e., the identity portion, ensure Q^l, 

Q^2 and Q^3 hold for ab &Vyj E. 

Proof. This lemma follows from direct substitution. In each case, there is only one nonzero element c^^ or 
''a(a+i) hand side of the equation from Q^l, Qg2 or Qg3 and substituting 1 for that coefficient gives 

the desired equality. □ 

It remains to define A', i.e. those coefficients c^-f, with ab G D and verify the corresponding equations in 

Qgl, Qc2, and Qj,3. For each column of A', Q^l, Qc2, and Qg3 yield a system of six scalar equations for 
the 2n variables {c^b^ijevvjE- Since we have many more variables than equations, there remains significant 
flexibility in the construction of a solution. In the upcoming sections, we will present such a solution where 
all but six of the coefficients in each column of A' are set to zero. The non-zero coefficients are chosen to be 
•-ol" ^''i "-ob' ^ol"^^'') c^}^ ^\ c^ab^ ^^'^ '^ab^^^ these have a natural correspondence to the geometry of the 
polygon and the edge ab; see Figure 4. 

We will show that the system of equations Q^l-Q^3 with this selection of non-zero coefficients for A' has 
an explicitly constructible solution. The solution is presented for special classes of polygons in Section 4 
and for generic convex polygons in Section 5. In each case, we prove a uniform bound on the size of the 
coefficients of A, a sufficient result to control as the following lemma shows. 
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Lemma 3.3. Let Q be a convex polygon satisfying Gl, G2. and G3, and let \i denote a set of barycentric 
coordinates satisfying the result of Theorem 2.3 (e.g. X^", A^'*'', A^'^'^'^, A^'*'*', or A^^'^'^. Suppose there exists 
M > 1 such that for all entries of A' , \c^f^\ < M. Then the functions are all bounded in H^, i.e. there 
exists a constant B > such that 

116,11^1(^2) (3-9) 
Proof. Since is defined by (3.1), Corollary 2.5 implies that there exists C > such that 

1 16,- 1 1 <C||A||. 

Since the space of linear transformations from E"("+i)/2 to M^" is finite-dimensional, all norms on A are 
equivalent. Thus, without loss of generality, we interpret ||A|| as the maximum absolute row sum norm, i.e. 

||A||:=max^|4|. (3.10) 

ab 

By the structure of A from (3.8) and the hypothesis, we have 

l|A||<^M 

□ 



4 Special Cases of the Serendipity Reduction 

Before showing that Q^l-Q^3 can be satisfied in a general setting, we study some simpler special cases in 
which symmetry reduces the number of equations that must be satisfied simultaneously. 

4.1 Unit Square 

We begin with the case where serendipity cilements were first examined, namely over meshes of squares. In 
recent work by Arnold and Awanou [3], the quadratic serendipity space on the unit square, denoted S2{I^), 
is defined as the span of eight monomials: 

^2(7^) ■.= span{l,x,y,x'^,xy,y'^,x^y,xy'^} (4.1) 

We will now show how our construction process recovers this same space of monomials. Denote vertices on 
[0,l]2by 

Vl = (0,0) V2 = (l,0) V3 = (l,l) V4 = (0,l) 
Vi2 = (1/2,0) V23 = (1,1/2) V34 = (1/2,1) Vi4 = (0,1/2) (4.2) 
Vis = V24 = (1/2, 1/2) 
The standard bilinear basis for the square is 

Ai = (l-x)(l-t/) A2 = a;(l-y) 

A4 = (1 - x)y A3 = xy 

Since the Aj have vanishing second derivatives and satisfy the definition of barycentric coordinates, they are 
in fact the harmonic coordinates A^*'' in this special case. Pairwise products give us the following 10 (not 
linearly independent) functions 

Mil = (1 - xf{l - yf /X12 = (1 - x)x{l - yf 

H22 = x'^{l - yf H23 = x'^{l - y)y 

M33 = /U34 = (1 - x)xy'^ 

H44 = {l-xfy'^ 1114 = {1 - x)^ (1 - y)y 

Hi3 = (1 - x)x{l - y)y H2i = (1 - x)x{l - y)y 



10 



For the special geometry of the square, = /i24, but this is not true for general quadrilaterals as we see 
in Section 4.3. The serendipity construction eliminates the functions fXis and /Z24 to give an 8-dimensional 
space. The basis reduction via the A matrix is given by 



" 61 " 




62 




63 




64 




62 




63 




64 




. 64 . 








1 

1 



































1 

1 
1 







Mil 


-1 







M22 





-1 




M33 


-1 







M44 





-1 




/il2 


1/2 


1/2 




M23 


1/2 


1/2 




M34 


1/2 


1/2 




Ml4 


1/2 


1/2 




Mis 






M24 



(4.3) 



It can be confirmed directly that (4.3) follows from the definitions of A given in the increasingly generic 
settings examined in Section 4.2, Section 4.3 and Section 5. The resulting functions are 



= {1 - x){l - y){l - X - y) 

C22 = x(l - y){x - y) 
£.33 = xy{-l + x + y) 
C44 = (1 - x)y{y - x) 



C12 = (1 - x)x{l - y) 
63 = x(\ - y)y 

C34 = (1 - x)xy 

iu = {l-x){l-y)y 



(4.4) 



Theorem 4.1. For the unit square, the basis functions {6i} defined in (4-4) satisfy Q^l-Q^3. 

Proof. A simple proof is to observe that the coefficients of the matrix in (4.3) satisfy Q^l-Q^3 and then apply 
Lemma 3.1. To illuminate the construction in this special case of common interest, we state some explicit 
calculations. The constant precision condition Q^l is verified by the calculation 

61 + 62 + Css + Ui + 262 + 263 + 2^4 + 2^4 = 1. 
The X component of the linear precision condition Q^2 is verified by the calculation 

(Vl)a;6l + (V2)x62 + (V3)x63 + (v4)xC44 + 

2(vi2)xCl2 + 2(v23)x63 + 2(V34)x64 + 2(vi4)x64 

= 62 + 63 + 2 • ^Ci2 + 2 • 1^3 + 2 • ^64 

= X. 

The verification for the y component is similar. The xy component of the quadratic precision condition Q^3 
is verified by 

(Vl)x(vi)j,^ll + (V2)x(v2)j^62 + (V3)x(v3)jy63 + (v4)x (v4)j;64 
+ [(Vl)x(v2)y + (V2)x(vi)y] ^12 + [(V2)x(v3)j, + {v 3) {v 2) y] ^23 
+ [(V3)x(v4)y + (V4)x(v3)j;] 64 + [(V4)x(vi)y + (vi)x(v4)j^] ^14 

= 63 + 63 + 64 = xy. 

The monomials x"^ and can be expressed as a linear combination of the ^ij similarly, via the formula given 
in Qj3. □ 

Corollary 4.2. The span of the functions defined by (4-4) is the standard serendipity space, i.e. 

span{^ij,^i(i+i)} =52(/^) 

Proof. Observe that x'^y = ^23 + £33 and xy"^ = ^33 + ^34. By the definition of S2{l'^) in (4.1) and the 
theorem, span {6i)6(i+i)} ^ S2{I^). Since both spaces are dimension eight, they are identical. □ 
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Figure 5: Notation for the construction for a regular polygon. 



4.2 Regulcir Polygons 

We now generalize our construction to any regular polygon with n vertices. Without loss of generality, this 
configuration can be described by two parameters 0<(T<6'<7r/2as shown in Figure 5. Note that the n 
vertices of the polygon are located at angles of the form ka where fc = 0, 1, . . . , n — 1. 

For two generic non-adjacent vertices and Vb, the coordinates of the six relevant vertices (recalling 
Figure 4) are: 



cos 6 
svaO 

cos 6 
— sva.6 



V5-1 = 



COS(t' — (7) 

sin{9 — a) 

cos{9 + a) 
- sin{9 + a) 



V6+1 = 



008(6* + cr) 
s'm(e + (j) 

cos{9 — a) 
— sm{9 — a) 



We seek to establish the existence of suitable constants c' 



ab' 



a,a+l a— l,a 



„bb 



6-1, b b,b+l 1- 1 

"'■ ' ^ab which 



„bb 
"'ab ' 



a— l,a 



preserve quadratic precision and to investigate the geometric conditions under which these constants become 
large. The symmetry of this configuration suggests that c'^^ = 
are reasonable requirements. For simplicity we will denote these constants by Cq := c^^, c_ := c^^^'", and 



b,b+l J o,o+l b-l,b 

Cnh > and c„: = c. 



c+ := 



a,a+l 
^ab 



Thus equation Q^l (which contains only six non-zero elements) reduces to: 

2co + 4c_ + 4c+ = 2. 



(4.5) 



Qg2 involves two equations, one of which is trivially satisfied in our symmetric configuration. Thus, the 
only restriction to maintain is 



2 cos 9co + 2 [cos 9 + cos{9 - a)] c_ + 2 [cos 9 + cos{9 + a)] c+ = 2 cos 9. 



(4.6) 



Qg3 gives three more requirements, one of which is again trivially satisfied. This gives two remaining 
restrictions: 



2 cos^ 9co +4 cos 9 cos{9 - o-)c_ + 4 cos 9 cos{9 + (t)c+ = 2 cos^ 9; 
2 sin^ 9cq + 4 sin 6* sin(6' - (t)c_ + 4 sin 9 sm{9 + ct)c+ = -2 sin^ ( 



(4.7) 
(4.8) 



Now we have four equations (4.5)-(4.8) and three unknowns cq, c_ and c+. Fortunately, equation (4.6) 
is a simple linear combination of (4.5) and(4.7); specifically (4.6) is times (4.5) plus 2^^^^ times (4.7). 
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With a little algebra, we can produce the system: 



1 2 2 

1 2(cosc7 + sincrtan^) 2(coscr — sin cr tan (?) 
1 2(cosa- — sinc7Cot6') 2(coscr + sincrcot ^) 

The solution of this system can be computed: 

(— 1 + cos ct) cot ^ + (1 + cos a) tan 6 







1 






1 


^+ . 




-1 



(4.9) 



Co 



(- 1 + cos cr) (cot ^ + tan 61) 



cos a — sin cr tan — 1 
2 (tan 6 + cot 6) sin a (cos a — 1)' 



1 — cos (T — sin (J tan 9 
2 (tan 6* + cot 0) sin cr (cos cr — 1) 



Although tan 6* (and thus the solution above) is not defined for 9 = it/2, the solution in this boundary 
case can be defined by the limiting value which always exists. We can now prove the following. 



Theorem 4.3. For any regular polygon, the basis functions {^ij} constructed using the coefficients 
= CO, = 4"+^ = c_, cr+^ = 4^'^ = c+ satisfy Q,l-Q,3. 

Proof. The construction above ensures that the solution satisfies Q^l, Q^.^, and Q^S. 



b 



□ 



The serendipity clement for regular polygons can be used for meshes consisting of only one regular 
polygon or a finite number of regular polygons. The former occurs only in meshes of triangles, squares 
and hexagons as these are the only regular polygons that can tile the plane. On the other hand, many 
tilings consisting of several regular polygons can be constructed using miiltiplc regular polygons. Examples 
include the snub square tiling (octagons and squares), the truncated hexagonal tiling (dodecahedra and 
triangles), the rhombitrihexagonal tiling (hexagons, squares, and triangles), and the truncated trihexagonal 
tiling (dodecagons, hexagons, and squares); see e.g. [6]. The construction process outlined above opens 
up the possibility of finite element methods applied over these types of mixed-geometry meshes, a mostly 
unexplored field. 



4.3 Generic Quadrilaterals 

Fix a convex quadrilateral with vertices Vi, V2, V3, and V4, ordered counterclockwise. We will describe 
how to set the coefficients of the submatrix A' in (3.8). It suffices to describe how to set the coefficients 
in the '13'-column of the matrix, i.e., those of the form cf'g. The '24'-column can be filled using the same 
construction after permuting the indices. Thus, without loss of generality, suppose that Vi := (—^,0) and 
V3 := {£, 0) so that V2 is below the x-axis and V4 is above the a;-axis, as shown in Figure 6. We have eight 
coefficients to set: 

11 22 33 44 12 23 34 , 14 

''131 ''13! ''13' ''13' ''13' ''13' ''13' auuCj^g. 

Using a subscript x or y to denote the corresponding component of a vertex, define the coefficients as follows. 

ctt:=0 (4.10) 



-13 

dl := ^ I""') ^ dt := , I""') , (4.11) 

(V4)y-(V2)j, (V2)s, - (V4)j, 

:= cll dt := dt (4.12) 

11 _ cli(V2)x + C?t(V4)c. ^ 33 _ c}i(V2)x+cf|(V4)x ^ 



, ^ ^13. , ^ (4.13) 

Note that there following the strategy shown in Figure 4, there are only six non-zero entries. For ease of 
notation in the rest of this section, we define the quantity 

, _ Cli(v2).+Cfi(v4). 

a.- ^ 

First we assert that the resulting basis does span all quadratic polynomials. 
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Figure 6: A generic convex quadrilateral, rotated so that one of its diagonals lies on the x-axis. Geometrically, 
cJI and cfl are selected to be coefficients of the convex combination of V2 and V4 that lies on the a;-axis. 



Theorem 4.4. For any quadrilateral, the basis functions constructed using the coefficients given in 

(4.10)-(4.13) satisfy Q^l-Q^3. 

Proof. Considering Lemmas 3.1 and 3.2, we only must verify Q^l-Qg3 in the cases when ab G D = {13, 24}. 

This will be verified directly by substituting (4.10)-(4.13) into the constraints Q^l-Q^S in the case ab = 13. 
As noted before, the ab = 24 case is identical, requiring only a permutation of indices. First note that 

ell + ell = -2 and cH - c?| = 2d. (4.14) 

For Q^l, the sum reduces to 

cll + cll + iicll + clt) = -2 + 4(1) = 2, 
as required. For Qg2, the a;-coordinate equation reduces to 

£{41 - cll) +2d£ = 

by (4.14) which is the desired inequality since we fixed (without loss of generality) Vab = (0,0). The y- 
coordinate equation reduces to 2(c}|(v2)j, + cf|(v4)y) = which holds by (4.11). Finally, a bit of algebra 
reduces the matrix equahty of Q^3 to only the equality ^^(cjg + cf|) = —2£'^ of its first entry (all other entries 
are zero), which holds by (4.14). □ 

Theorem 4.5. Over all convex quadrilaterals, ||A|| is uniformly bounded. 

Proof. By Lemma 3.3, it suffices to bound Ic^^l uniformly. First observe that the convex combination of the 
vertices V2 and V4 using coefficients cJI and cff produces a point lying on the a;-axis, i.e., 

1 =cll + cll, and (4.15) 

0=c\l{V2)y + 4i{V4)y. (4.16) 

Since (v2)j; > and (v4)y < 0, (4.12) implies that cll,clj G (0,1)- By (4.12), it also follows that cf|,c}| € 
(0,1). 

For c 1 and cfl, note that the quantity d£ is the a;- intercept of the line segment connecting V2 and V4. Thus 
d£ G [-£, £] by convexity So d e [-1,1] and thus (4.13) implies \c\l\ = \d-l\ < 2 and |cf|| = \-d-l\ < 2. □ 
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Figure 7: Generic convex polygon, rotated so that Vo = (—-^,0) and = (£,0). The x-intercept of the 
line between Va-i and Va+i is defined to be —da£ and the a;-intercept of the line between v^-i and Vb+i is 
defined to be dbi. 



5 Proof of the Serendipity Reduction on Generic Convex Polygons 

We now define the sub-matrix A' from (3.8) in the case of a generic polygon. Pick a column of A', i.e., fix 
ab € D. The coefficients cj'^ are constrained by a total of six equations Q^l, Qg2, and As before (recall 
Figure 4), six non-zero coefficients will be selected in each column to satisfy these constraints. Specifically, 

c^b := 0, for i {a, b} and 0^^^+^^ = 0, for i ^ {a - l,a,b - 1, b}, (5.1) 
leaving only the following six coefficients to be determined: 

,66 (a-l)a a(a+l) (6-1)6 „^ , 6(6+1) 



aa Ob AU-i-ja aya^i) \u-i,u , 



ab 



For the remainder of this section, we will omit the subscript ab to ease the notation. Writing out Q^l-Q^B 
for this fixed ab pair, we have six equations with six unknowns: 

^aa _,. ^bb _,. 2c(«-i)" + 2c''(°+^) + 2c(*'-i)*' + 2c^(*+i) = 2; 

c'^Vaa + 2c(»-l)V(,_i)„ + 2C«("+1)V,(,+ 1) + 

€"^66 + 2c(''-l)V(6_i)6 + 2c''(''+l)v6(6+i) = 2v„6; 

c-v.vr + c(«-i)»(v„_ivr + v„vri) + c"("+i)(v„vr+i + v„+ivr)+ 

c'-^fevf + c(''-l)''(v6-ivf + V^V^i) + C''(''+1)(V6V^+1 + V6+lVr) = V.V^ + V6V^ 

Assume without loss of generality that = {—£, 0) and V6 = {£, 0) with £ < 1/2 (since Q has diameter 
1). We introduce the terms da and df, defined by 

da := K-i)^K+i), - (v,+l).(Va-l), . 1 ^5 2) 

(Va-l)y - (Va+l)y ^ 
(V6+l)a^(v6-l)^ - (v6-l)a;(v6+l)y 1 



15 



These terms have a concrete geometrical interpretation as shown in Figure 7: —da£ is the x-intercept of 
the hne between Va-i and Vo+i, while db£ is the a;-intercept of the line between Vb_i and Vb+i. Thus, 
by the convexity assumption, da,rf(, € [—1,1]- Additionally, —da < db with equality only in the case of a 
quadrilateral which was dealt with previously. For ease of notation and subsequent explanation, we also 
define 

First we choose c^""^)" and c°'(°'~^^^ as the solution to the following system of equations: 

^(a-l)a ^ ^a{a+l) ^ (5 5) 

c(»-l)V„_i + c"("+^)va+i = sda^Ta. (5.6) 

There are a total of three equations since (5.6) equates vectors, but it can be verified directly that this system 
of equations is only rank two. Moreover, any two of the equations from (5.5) and (5.6) suffice to give the 
same unique solution for c^"--^)°- and c"("+^'. 

Similarly, we select c^**"^)^ and c*'^''"^) as the solution to the system: 



c 



{b-i)b. 



Finally, we assign c"" and c"" by 



V6-1 + c''(''+iV6+i = sdbVb- (5.8) 



^'^'^ and (5.9) 



2-{da + db) 

and claim that this set of coefficients leads to a basis with quadratic precision. 

Theorem 5.1. For any convex polygon, the basis functions {^ij} constructed using the coefficients defined 

by (5.5)- (5. 10) satisfy Q^l-Q^S. 

Proof. Based on Lemmas 3.1 and 3.2, it only remains to verify that Q^l, Qc2, and Qg3 hold when ab £ D. 
Observe that c"" and (^'' satisfy the following equations: 

+ c'* + 4s = 2; (5.11) 
c«« - c*-^ + - 4) = 0; (5.12) 

^aa ^^bb ^2s{da + db) = -2. (5.13) 

First, note that Q^l follows immediately from (5.5), (5.7) and (5.11). 

The linear precision conditions (Qc2) are just a matter of algebra. Equations (5.5)-(5.8) yield 

C""V„„ + C^^'Wbb + 2c(«-l)°V(„_i)„ + 2c«(«+l)v„(„+i) + 2c(''-l)''V(6_i)6 + 2c*(*'+l)vb(6+i) 

= (C«" + c(«-^)" + c"(«+l))Va + {C^^ + + C^^^+^^)Vb + sdaVa + sdbVb 

= {C"" + S + sda)Va + {c^^ + S + sdb)Vb. 

Substituting the fixed coordinates of Va = {—i, 0) and = {£, 0) reduces this expression to the vector 

(-c"« -s-sda + c'>'' + s + sdb)e 


Finally, we address Qc3. Factoring the left side gives. 
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Figure 8: Notation used in proof of Theorem 5.3. 



a(a+l) T 



a+l) 



= c""VaV^ + C^^b^b + sdaVa^l + sdbVbVb + Va(srfaVa ) + Vb{sdbVb) 
= (C»» + 2s4)VaV^ + (C*-^ + 2s4)v6Vr 

Again substituting the coordinates of and v?,, we obtain the matrix 



2sda + C 




bh 



The right side of Q^S is 



v„v. 



T 







Hence the only equation that must be satisfied is exactly (5.13). 



□ 



Remark 5.2. We note that s was specifically chosen so that (5.11)-(5.13) would hold. The case s = 1 happens 
when da = — dfe, i.e. only for the quadrilateral. 



Theorem 5.3. Given a convex polygon satisfying Gl, G2 and G3, 



is uniformly hounded. 



Proof. By Lemma 3.3, it suffices to show a uniform bound on the six coeflacients defined by equations 
(5.5)-(5.10). First we prove a uniform bound on da and db given G1-G3. 

Wc fix some notation as shown in Figure 8. Let C(va,(i*) be the circle of radius (from G2) around 
Va. Let p~ := {Px,Py) and p+ := {pt^Py) be the points on C(vo,d*) where the line segments to from 
Va_i and Va+i, respectively, intersect. The chord on C{va,d^) between p~ and p+ intersects the a;-axis at 
Xp := (a;p,0). By convexity, {va)x < Xp. 
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To bound Xp — (va)^; below, note that the triangle v^p p+ with angle /?„ at Vq is isosceles. Thus, the 
triangle VaP^Xp has angle Zv^p^Xp = (tt — /3a)/2, as shown at the right ol Figure 8. The distance to the 

nearest point on the line segment between p" and p+ is d* sin ^ ^^2" ) ■ Based on G3, e« > is defined to be 

Xp - {^a)x > sin (—^^^ > ^* (^ ^ 2 ) ^" ^* ^ ^" '^■^'^ 

Since —dj. < 1 is the .x-intercept ol the line between Va_i and v^+i, we have Xp < —dai. Then we rewrite 
(va)a; — in the geometrically suggestive form 

{Xp - {Wa)x) + {-dj - Xp) + (0 + dj) = i. 

Since —daf. — Xp > 0, we have Xp — {va)x + da( < (■ Using (5.14), this becomes da£ < ^ — e*. Recall from 
Figure 7 and previous discussion that da, df, € [—1, 1] and —da < db- By symmetry, di,£ < £ — e* and hence 
da + db< 2^-2e* < l-2e*. 

We use the definition of c"" from (5.9), the derived bounds on da and db, and the fact that £ < 1/2 to 
conclude that 

,„aa, . |2 + 2rf„| ^ 2 + 2(1 - je./e)) ^ 4 - 46, ^ 



' ' 1 + 2e, 1 + 2e* " 1 + 2e, 

Similarly, |c^''| < < 4. For the remaining coefficients, observe that the definition of s in (5.4) implies 

that < s < 2/(1 + 2e,). Equation (5.5) and the y-component of equation (5.6) ensure that c*^"~^)"/s and 
^a{a+i)^g are the coefiicients of a convex combination of Va-i and v^+i. Thus c("^i)'', 0^°+^^" e (0, s) and s 
serves as an upper bound on the norms of each coefficient. Likewise, |c(*'~^)*'|, |c''(^"'"^)| < s. Therefore, 

/4^4e^ 2 

max , , 1 

Vl + 2e*' l + 2e*' 

is a uniform bound on all the coefficients of A. □ 



6 Converting Serendipity Elements to Lagrange-like Elements 

The 2n basis functions constructed thus far naturally correspond to vertices and edges of the polygon, but 
the functions associated to midpoints are not Lagrange-like. This is due to the fact that functions of the form 
may not evaluate to 1 at Vj(j+i) or may not evaluate to at Vj(j_,_i), even though the set of {^ij} 
satisfies the partition of unity property Qjl. To fix this, we apply a simple bounded linear transformation 
given by the matrix B defined below. 

To motivate our approach, we first consider a simpler setting: polynomial bases over the unit segment 
[0,1] C M. The baryccntric functions on this domain are Xo{x) = 1 — x, and Ai(a;) = x. Taking pairwise 
products, we get the quadratic basis fJ.oo{x) '■= (Ao(a;))^ = (1 — x)^, ijqi{x) := Xo{x)Xi{x) = (1 — x)x, and 
^11 (x) :— (Ai(x))^ = x^, shown on the left of Figure 9. This basis is not Lagrange-like since ^01 (1/2) 1 
and /ioo(l/2), /iii(l/2) 7^ 0. The quadratic Lagrange basis is given by ipoo{x) '■— 2(1 — x) (| — x), ipoi{x) '■= 
4(1 — x)x, and t/'ii(x) := 2 (x — |) x, shown on the right of Figure 9. These two bases are related by the 
linear transformation Bid: 





ipoo 




"1 





-f 








i>ii 







1 


-1 




Mil 




.V'Ol. 










4 




Moi 



^Bid[m, 



(6.1) 



This procedure generalizes to the case of converting the 2D serendipity basis {^ij} to a Lagrange like 
basis {i>ij}- Define 

V'ii := ^ii - - ?i-i,i and tpi,i+i = 4 
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Figure 9: A comparison of the product barycentric basis (left) with the standard Lagrange basis (right) for 
quadratic polynomials in one dimension. 

Using our conventions for basis ordering and index notation, the transformation matrix B taking [^ij] to [tpij] 
has the structure 



■023 



'nl 



-1 -1 



^12 
inl. 



The following proposition says that the functions {tpij} defined by the above transformation are Lagrange- 
likc. 

Proposition 6.1. For all i,j n}, ipuivj) = S^, tpuivjj+i) = 0, V'i(i+i)(vj) = 0, and V'i(i+i)(vjj+i) = 

Proof. We show the last claim first. By the definitions of B and A, we have 



a<b 



Since Xj is picccwisc linear on the boundary of the polygon, Aa(vjj-|_i) = 1/2 if a G {j, j + 1} and zero 
otherwise. Accordingly, Haai'^jj+i) = 1/4 if a € {j,j + 1} and zero otherwise, while ^ab{'^j,j+i) = 1/4 if 
{a, b} = + 1} and zero otherwise. 



since the identity structure of A as given in (3.8) implies that Cj^j^^^ = c^j'+i)^(j_|_i) = and that Cj^^^i) = ^ij- 
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Next, observe that fiabi^j) = ^a{^j)^b{^j) — 1 ii a = b = j and otherwise. Hence, any term of the 
form c*^/iab(vj) for a 7^ 6 is necessarily zero. Therefore, by a similar expansion, ^j(j_|_i)(vj) = Cj^j~^^^ = 0, 
proving the penultimate claim. 

For the first two claims, similar analysis yields 

^iii^j) = 6i(Vj) - 6(i+l)(Vj) - ^(^-l)i(Vj) 
— J'^ .-I _ J(»+l) . 1 _ . 1 

again by the identity structure of A. Finally, by similar analysis, wc have that 

= - Ci(i+i)(vj,j+i) - ^(i-i)i(vjj+i) 

= {Sij + SiQ+i))- - -Sij - -di(j+i) = 0, 

completing the proof. □ 
In closing, note that ||B|| is uniformly bounded since its entries all lie in {—1,0, 1,4}. 



7 Applications and Extensions 

Our quadratic serendipity element construction has a number of uses in modern finite element application 
contexts. First, the construction for quadrilaterals given in Section 4.3 allows for quadratic order methods on 
arbitrary quadrilateral meshes with only eight basis functions per element instead of the nine used in a bilinear 
map of the biquadratic tensor product basis on a square. In particular, we show that our approach maintains 
quadratic convergence on a mesh of convex quadrilaterals known to result in only linear convergence when 
traditional serendipity elements are mapped non-affinely [4]. 

We solve Poisson's equation on a square domain composed of trapezoidal elements as shown in Fig- 
ure 10 (left). Boundary conditions are prescribed according to the solution u{x,y) = sin(a;)e^; we use our 
construction from Section 4.3 starting with mean value coordinates {A^^'^^}. Mean value coordinates were 
selected based on a few advantages they have over other types: they are easy to compute based on an explicit 
formula and the coordinate gradients do not degrade based on large interior angles [27]. For this particular 




n 



n 



4 





u — -ah 




\\V{u- 


Uh)\\i^1 


n 


error 


rate 


error 


rate 


2 


2.34e-3 




2.22e-2 




4 


3.03e-4 


2.95 


6.10e-3 


1.87 


8 


3.87C-5 


2.97 


1.59e-3 


1.94 


16 


4.88C-6 


2.99 


4.04e-4 


1.97 


32 


6.13C-7 


3.00 


1.02e-4 


1.99 


64 


7.67e-8 


3.00 


2.56e-5 


1.99 


128 


9.59e-9 


3.00 


6.40e-6 


2.00 


256 


1.20e-9 


3.00 


1.64e-6 


1.96 



Figure 10: Trapezoidal meshes (left) fail to produce quadratic convergence with traditional serendipity 
elements; see [4]. Since our construction begins with affinely-invariant generalized barycentric functions, the 
expected quadratic convergence rate can be recovered (right). The results shown were generated using the 
basis {t/'ij} resulting from the selection of the mean value coordinates as the initial barycentric functions. 
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example, where no interior angles asymptotically approach 180°, Waclispress coordinates give very similar 
results. As shown in Figure 10 (right), the expected convergence rates from our theoretical analysis are 
observed, namely, cubic in the i^-norm and quadratic in the i7^-norm. 

An additional application of our method is to adaptive finite elements, such as the one shown in Figure 11. 
This is possible since the result of Theorem 5.3 still holds if G3 fails to hold only on a set of consecutive 
vertices of the polygon. This weakened condition suffices since consecutive large angles in the polygon do 
not cause the coefficients cj'^ to blow up. For instance, consider the degenerate pentagon shown in Figure 11 
which satisfies this weaker condition but not G3. Examining the potentially problematic coefficients C25, 
observe that the lines through vi, V4 and Vi, V3 both intersect the midpoint of the line through V2, V5 
(which happens to be Vi). In the computation of the C25 coefficients, the associated values d2 and are 
both zero and hence s = 1 (recall Figure 7 and formula (5.4)). Since s is bounded away from 00, the analysis 
from the proof of Theorem 5.3 holds as stated for these coefficients and hence for the entire element. A more 
detailed analysis of such large-angle elements is an open question for future study. 




Figure 11: Theorem 5.3 can be generalized to allow certain types of geometries that do not satisfy G3. The 
degenerate pentagon (left), widely used in adaptive finite element methods for quadrilateral meshes, satisfies 
Gl and G2, but only satisfies G3 for four of its vertices 

still hold on this geometry, resulting in the Lagrange-like quadratic element (right). 



The bounds on the coefficients c^^^ from Section 5 



Nevertheless, the geometric hypotheses of Theorem 5.3 cannot be relaxed entirely. Arbitrarily large 
non-consecutive large angles as well as very short edges, can cause a blowup in the coefficients used in 
the construction of A, as shown in Figure 12. In the left figure, as edges Va_iVa and Vb-iV;, approach 
length zero, da and dh both approach one, meaning s (in the construction of Section 5) approaches 00. 



In this case, the coefficients 



-'ij 



and ^•'^ grow larger without bound, thereby violating the result of 



Theorem 5.3. In the right figure, as the overall shape approaches a square, da and db again approach one 



so that s again approaches 00. In this case, all the coefficients c. 



(a— l)a 

ij 



a(a+l) 



and c^j^"*"^^ all grow 



without bound. Nevertheless, if these types of extreme geometries arc required, it may be possible to devise 
alternative definitions of the coefficients satisfying Q^l-Q^S with controlled norm estimates since the set 
of restrictions Qgl-Qg3 does not have full rank. Note that this flexibility has lead to multiple constructions of 
the traditional serendipity square [24, 22]. Cursory numerical experimentation suggests that some bounded 
construction exists even in the degenerate situation. 

The computational cost of our method is an important consideration to application contexts. A typical 
finite element method using our approach would involve the following steps: (1) selecting A; coordinates and 
implementing the corresponding tpij basis functions, (2) defining a quadrature rule for each affine- equivalent 
class of shapes appearing in the domain mesh, (3) assembling a matrix L representing the discrete version 
of linear operator, and (4) solving a linear system of the form L?i = /. The quadrature step may incur 
some computational effort, however, if only a few shape templates are needed, this is a one-time fixed 
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Figure 12: The hypotheses of Theorem 5.3 cannot be relaxed entirely as demonstrated by these shapes. If 
G2 does not hold, arbitrarily small edges can cause a blowup in the coefficients c^^^ (left). If G3 does not 
hold, non-consecutive angles approaching tt can cause a similar blowup. 



pre-preprocessing cost. In the trapezoidal mesh example from Figure 10, for instance, we only needed 
one quadrature rule as all domain shapes were affinity equivalent. Assembling the matrix L may also be 
expensive as the entries involve integrals of products of gradients of ^pij functions. Again, however, this cost 
is incurred only once per affine-equivalent domain shape and thus can be reasonable to allow, depending on 
the application context. 

The computational advantage to our approach comes in the final linear solve. The size of the matrix L is 
proportional to the number of edges in the mesh, matching the size of the corresponding matrix for quadratic 
Lagrange elements on triangles or quadratic serendipity elements on squares. If the pairwisc products ^ab 
were used instead of the tpij functions, the size of L would be proportional to the square of the number of 
edges in the mesh, a substantial difference. 

Finally, we note that although this construction is specific to quadratic elements, the approach seems 
adaptable, with some effort, to the constriiction of cubic and higher order serendipity elements on generic 
convex polygons. As a larger linear system must be satisfied, stating an explicit solution becomes complex. 
Further research along these lines should probably assert the existence of a uniformly bounded solution 
without specifying the construction. In practice, a least squares solver could be used to construct such a 
basis numerically. 
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